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Abstract — Much research has been devoted to the problem of 
restoring Poissonian images, namely for medical and astronom- 
ical applications. However, the restoration of these images using 
state-of-the-art regularizers (such as those based on multiscale 
\ representations or total variation) is still an active research area, 
since the associated optimization problems are quite challenging. 
In this paper, we propose an approach to deconvolving Poissonian 
images, which is based on an alternating direction optimization 
method. The standard regularization (or maximum a posteriori) 
restoration criterion, which combines the Poisson log-likelihood 
with a (non-smooth) convex regularizer (log-prior), leads to hard 
optimization problems: the log-likelihood is non-quadratic and 
non-separable, the regularizer is non-smooth, and there is a 
non-negativity constraint. Using standard convex analysis tools, 
we present sufficient conditions for existence and uniqueness of 
solutions of these optimization problems, for several types of 
\ regularizers: total-variation, frame-based analysis, and frame- 
, based synthesis. We attack these problems with an instance of 
the alternating direction method of multipliers (ADMM), which 
belongs to the family of augmented Lagrangian algorithms. We 
study sufficient conditions for convergence and show that these 
are satisfied, either under total-variation or frame-based (analysis 
and synthesis) regularization. The resulting algorithms are shown 
to outperform alternative state-of-the-art methods, both in terms 
of speed and restoration accuracy. 

Index Terms — Image restoration, image deconvolution, Pois- 
son images, convex optimization, alternating direction methods, 
augmented Lagrangian. 



I. Introduction 

A large fraction of (if not all) the work on image denoising, 
restoration, and reconstruction has been devoted to developing 
regularizers (priors, from a Bayesian point of view) to deal 
with the presence of noise and/or the ill-conditioned or ill- 
posed nature of the observation operator, and to devising 
efficient algorithms to solve the resulting optimization prob- 
lems. Much of that work assumes linearity of the observation 
operator {e.g., the convolution with some point spread func- 
tion, the acquisition of tomographic projections, or simply an 
identity in the case of denoising) and the presence of additive 
Gaussian noise. For this classical scenario, recent state-of- 
the-art methods adopt non-smooth convex regularizers, such 
as total-variation or the li norm of frame coefficients; the 
resulting optimization problems are convex but non-smooth, of 
very high dimensionality, and have stimulated a considerable 
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amount of research on special purpose algorithms (see fllj El, 
|9|, 1 18 1, 1 40 1, 1 45 1 and the many references therein). 

The algorithms developed for linear operators and Gaussian 
noise cannot be directly applied to other observation models, 
such as the Poisson case considered in this paper Poissonian 
image models are well studied and highly relevant in fields 
such as astronomical |38|, biomedical mil, |[16l, EOl, ||33l . 
||4T|. p3l, and photographic imaging |fT9l . A very recent 
and comprehensive overview of deconvolution methods for 
Poissonian images can be found in fT2| (where a state-of- 
the-art algorithm is also introduced); we refer the reader to 
that publication for more references on this topic. 

The standard criterion for deconvolution of Poissonian im- 
ages consist of a convex constrained optimization problem: 
the objective function includes the so-called data term, which 
is convex and smooth, but not quadratic, plus a convex non- 
smooth regularizer (the log-likelihood and log-prior, from 
a Bayesian inference perspective), and a constraint forcing 
the solution to be non-negative. Although the problem is 
convex, its very high dimensionality (when dealing with 
images) usually rules out the direct application of off-the-shelf 
optimization algorithms. 

Furthermore, the Poisson log-likelihood, which is non- 
quadratic and non-separable (except in the pure denoising 
case) raises several difficulties to the current state-of-the- 
art image deconvolution algorithms. More specifically, the 
Poisson log-likelihood does not have a Lipschitz-continuous 
gradient, a sufficient condition for the applicability (with guar- 
anteed convergence) of algorithms of the forward-backward 
splitting (FBS) family ID, Oil, ES). If, nevertheless, an FBS 
method is applied, it is well known to be slow, specially if the 
observation operator is severely ill-conditioned, a fact which 
has stimulated recent research on faster methods HI, |[2], B31 ; 
these faster algorithms also require the log-likelihood to have 
a Lipschitz-continuous gradient, which is not the case with 
Poissonian observations. 

In this paper, we propose a new approach to tackle the 
optimization problem referred to in the previous paragraph. 
The proposed algorithm is based on an instance of the alter- 
nating direction method of multipliers (ADMM) II13I . II20L 
||2l|, which belongs to the family of augmented Lagrangian 
methods |[29l . For this reason, we call it PIDAL (Poisson 
image deconvolution by augmented Lagrangian). Although the 
proposed approach is related to the recent split-Bregman (SB) 
technique Il22l . our splitting strategy and resulting algorithm 
are quite different from the one in 122] (which, moreover, 
is not adequate for Poissonian image models). Finally, we 
mention that this paper is an extension of our much shorter 
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and preliminary work IfTTl . 

In recent work, Douglas-Rachford splitting (DRS) methods 
were proposed to attack problems in which log-likelihood the 
does not have a Lipschitz-continuous gradient lIS). In fact, the 
ADMM is closely related to DRS methods (13], (15], so the 
method proposed in this paper can also be interpreted from a 
DRS viewpoint. 

In this paper, we will consider three types of regulariza- 
tion: total variation |4|, [31] and both frame-based analysis 
and frame-based synthesis formulations llT4ll . In Section [III 
after presenting these three formulations, we derive sufficient 
conditions for existence and uniqueness of solutions of the 
corresponding optimization problems. The ADMM framework 
is reviewed in Section Hill where we also introduce the partic- 
ular variant that is suitable for a linear combination of several 
convex functions, which is the form of the objective function 
in hand. In Sections |IV] and |V] we instantiate the proposed 
variant of ADMM to the three types of regularizers considered 
and provide sufficient conditions for convergence. Finally, 
the effectiveness of the resulting algorithm is illustrated in 
comparison with current state-of-the-art alternatives [11 1, l,12l . 
jm, ll37l . Il39l . via a set of experiments reported in Section IVll 

II. Problem Formulation 

In this section, we begin by reviewing the derivation of 
the standard log-likelihood resulting from assuming that the 
observations are Poisson distributed with a mean intensity 
linearly related with the underlying image to be estimated. 
Then, we present three different regularization/Bayesian crite- 
ria, using synthesis and analysis formulations llT4l . and study 
existence/uniqueness of the corresponding solutions. 

A. The Linear/Poisson Observation Model 

Let y = (yi, ym) £ N™ denote an m-vector of observed 
counts (No = N U {0}), assumed to be a sample of a random 
vector Y — (Yi, Ym) of rn independent Poisson variables, 
with probability distribution 



P[Y = y|A] - n 



xr 



(1) 



where A = (Ai,...,Am) G K"' (IR+ denotes the non-negative 
reals) is the underlying mean (intensity) vector, assumed to be 
a linear observation of an unknown image x, i.e.. 



A = K: 



(2) 



where K the observation operator, which in our finite dimen- 
sional setting is simply a matrix K G R™^". This matrix 
may model a convolution or some other linear observation 
mechanism, such as emission tomography. So that the under- 
lying unknown x can also have the meaning of intensity, it is 
commonly assumed that x € R" . It is usually further assumed 
that all the elements of K are non-negative ifTTj . |[T6l . Il38l . 
When dealing with images, we adopt the usual vector notation 
obtained by stacking the pixels into a vector, in lexicographic 
order. 



Combining ([U and ^ and taking logarithms leads to the 
negative log-likelihood function ifTTI . ll38l . 



logP[Y = y|x] 



E 

1=1 



(Kx),-yaog((Kx) 



- £(Kx) 

= (^°K)(x), 



(3) 
(4) 



where (v)^ (or Vi) denotes the i-th component of some vector 
V and £ : R™ R = R U {-oo, +00} is the negative log- 
likelihood function for the case K = I, that is 



(5) 



Dealing with the particular case Zi — requires some care, 
because of the presence of the logarithm. Seen as function of 
z to be used in a minimization problem, it is convenient to 
write the negative log-likelihood function as 



(6) 



4=1 



where C = log(yi!) is a finite (recall that 0! = 1) irrelevant 
(independent of z) constant and ^ : R x Nq R is defined as 

^{z,y) ^ z + LM^{z) -y \og{z+), (7) 

where ls is the indicator function of set S, 

^ zeS 

+00 4= z ^ S, 



z+ = max{0,z}, log(O) = —00, and Olog(O) = 0. 

The following two propositions characterize ^ as a function 
of its first argument, as well as C and £ o K, in terms of the 
key concepts of convex analysis (see Appendix A). 

Proposition 1: For any y G Nq the function ^{■■,y) : R — > R 
is proper, lower semi-continuous (Isc), coercive, and convex. If 
y > 0, then ^(•. y) is also strictly convex. 

Proof: For y = 0, ^(z^O) = z + iR+(z), thus ^(-,0) 
is the sum of the identity function with tR^, which are 
both proper, Isc, coercive, and convex. For any y > 0, 
£,{z,y) = z + tR^(z) - ylog(z+); since ylog((-)+) is also 
proper, Isc, coercive, and convex, so is y). Finally, if y > 0, 
ylog((-)+) is strictly convex (see the definition in Appendix 
A), thus ^{■■,y) is also strictly convex. ■ 

Proposition 2: Function C is proper, Isc, coercive, and con- 
vex. If yi 7^ 0, for i = 1, ...,m, £ is also strictly convex. 
Function C oK. is proper, Isc, and convex. Function C o K. is 
coercive if K is injective. Function £ o K is strictly convex if 
K is injective and yt 7^ 0, for i = 1, ...,m. 

Proof: Function C is the sum of proper, Isc, coercive, 
convex functions. If yi ^ 0, for i — 1. m, the functions in 
the sum are also strictly convex, thus C is also strictly convex. 
Function £ o K is the composition of a proper, Isc, convex 
functions with a linear function, thus it is proper, Isc, and 
convex. If K is injective, its null set is the zero vector, thus 
lim||x||_^+oo l|Kx|| = +00, thus £ o K is coercive. Finally, 
if K is injective and y^ 7^ 0, for i = l,...,m, £ is strictly 
convex, thus so is £ o K. ■ 
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B. Regularization Criteria: Analysis and Synthesis Formula- 
tions 

Under a maximum a posteriori (MAP) or regularization cri- 
terion, the image estimate is obtained by solving a variational 
problem: minimizing an objective function, which includes the 
log-likelihood term plus a regularizer ifTTl . Il26l . Il38l . under 
a positivity constraint. We will now describe three possible 
ways of building such an objective function. 

1) Total Variation Regularization: A standard choice for 
regularization of digital image restoration/reconstruction prob- 
lems is the isotropic discrete TV regularizer [4], 

" I 

Tv(x) ^ V (^'^)' + i^i^y^ (8) 

s=l 

where (A^x and A^x) denote the horizontal and vertical first 
order differences at pixel s, respectively. This regularizer is a 
discrete version of the TV regularizer proposed in ifsTl . The 
resulting optimization problem is 

min L'^^(x) (9) 

X 

with 

iTV(x) = CiK x) + T TV(x) + (x), (10) 

where r e M+ is the regularization parameter and the role 
of iRn, the indicator of the first orthant, is to impose the 
non-negativity constraint on the estimate. The next proposition 
concerns the existence and uniqueness of minimizers of L^^ . 

Proposition 3: Consider the function L^^ defined in ( fTol i. 

a) L^^ is proper, Isc, and convex. 

b) If the intersection of the space of constant images 
{x = a(l, 1, 1), a e M} (which is the null space 
of TV) with the null space of K is just the zero 
vector, then L^^ is coercive, and (|9]l has a solution. 

c) If K is injective, then (|9]l has a solution. 

d) If K G and at least one element of K is 
strictly positive, then (|9]l has a solution. 

e) If K is injective and yi ^ 0, for i = 1, m, then 
L^^ is coercive and strictly convex thus there is a 
unique solution. 

Proof: 

a) The functions ir" , TV, and £ o K (Proposition |2]l 
are proper, Isc, and convex, thus so is their sum. 

b) Similar to 13. 

c) If K is injective, its null space is just the zero vector, 
thus £ o K and L^^ are coercive. 

d) If aU the elements of K are non-negative and at least 
one is positive, then the constant vector (1, 1, 1) 
doesn't belong to the null space of K and the result 
follows from (b). 

e) If K is injective and ^ 0, for i — l,...,m, £oK 
is strictly convex (Proposition |2]l, thus so is L^^ and 
its minimizer is unique. 



2) Frame Analysis Regularization: The use of a regu- 
larizer which is a direct function of the unknown image 
(as in (l9]l-(fT0ll) corresponds to a so-called analysis-based 
prior/regularizer llT4l . Another well-known type of analysis- 
based regularization penalizes the norm (typically £i) of the 
representation coefficients of x on some wavelet basis or 
tight frame, given by Px, where P is the analysis operator 
associated with the frame |28|. This approach leads to the 
following optimization problem: 

minL'''^(x), (11) 

X 

where FA stands for frame analysis and 

LP^(x) =/:(Kx) + r||Px||i+tMn(x); (12) 

as above, r is the regularization parameter and ir^ imposes the 
non-negativity constraint on the estimate. The next proposition 
addresses the existence and uniqueness of minimizers of L^^. 

Proposition 4: Consider the function L^^ defined in (fT2] l. 

a) L^'^ is proper, Isc, convex, and coercive, thus has a 
minimizer 

b) If K is injective and yi ^ 0, for i — 1, ...,m, then 
L^^ is strictly convex with a unique minimizer. 

Proof: 

a) The functions iR^, || • ||i oP, and £oK (Proposition 
|2]i are proper, Isc, and convex, thus so is their sum. 
Furthermore, since P is the analysis operator of a 
tight frame, its null space is simply the zero vector, 
thus II • 111 o P is coercive. 

b) If K is injective and yi ^ 0, for i = 1, ...,m, >CoK 
is strictly convex (Proposition |2]), thus so is L^^ and 
its minimizer is unique. 

■ 

3) Frame Synthesis Regularization: Finally, another well- 
known class of approaches is known as synthesis-based lfT4ll . 
Here, the unknown image is represented on a frame (e.g., 
of wavelets, curvelets, or other multi-scale system) and then 
the coefficients of this representation are estimated from the 
observed data, under some regularizer With W G K"'*'' 
denoting the synthesis matrix of the frame, the image is 
written as x = W s, where s is the vector of representation 
coefficients, and the resulting optimization problem is 

min L^^(s) (13) 

S 

where FS stands for frame synthesis and 

iPS(s) =/:(KWs)+t||s||i + 6r^(Ws). (14) 

Naturally, the indicator function tR^ forcing the image esti- 
mate to be non-negative is applied to the image W s and not 
its coefficients. The next proposition addresses the existence 
and uniqueness of minimizers of L^^. 

Proposition 5: Consider the function L^^ defined in ( IT4b . 

a) L^^ is proper, Isc, convex, and coercive, thus has a 
minimizer. 

b) If KW is injective and yi ^ 0, for i = l,...,m, 
then LF^ is strictly convex with a unique minimizer. 
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Proof: 

a) The functions if^^ o W, || • ||i, and C o KW 
(Proposition 121) are proper, Isc, and convex, thus so 
is their sum. Furthermore, since || • ||i is coercive, 
L^^ is coercive. 

b) Same proof as that of Proposition |4](b). 



III. The Alternating Direction Method of 
Multipliers (ADMM) 

A. The Standard ADMM 

The key tool in this paper is the alternating direction 
method of multipliers (ADMM) (|T3l, lEO), (IT]. Consider an 
unconstrained problem of the form 



min /i(z)+/2(Gz) 



(15) 



where fi : ^ M, f2 : W ^ M, and G e W^'K The 
ADMM for this problem is defined in Fig. [T] 



Algorithm ADMM 

1. Set k = Q, choose /i > 0, uo, and do. 

2. repeat 

3. Zfe+i e argminz /i(z) + f||Gz - Ufe - dfc||2 

4. Ufc+i G argminu/2(u) + |||Gzfe+i - u - dfc||i 

5. dfe+i dfc - (Gzfc+i - Ufc+i) 

6. fc <~ + 1 

7. until stopping criterion is satisfied. 



Fig. 1. The ADMM algorithm. 

For later reference, we now recall a theorem by Eckstein 
and Bertsekas in which convergence of (a generalized version 
of) ADMM is shown. 

Theorem 1 (Eckstein-Bertsekas, l[13]I ): Consider problem 
dTsl ). where G S MP^'^ has full column rank and 
/i : R'^ ^ R and f2 : Rp ^ M are closed, proper, 
convex functions. Consider arbitrary /i > and uq, dg G W. 
Let {rjk > 0, fc = 0, 1, ...} and {pk > 0, fc = 0, 1, ...} be two 
sequences such that 



r]k < oo and 

k=0 



^Pfe < oo. 

k=Q 



Consider three sequences {zk £ M"*, k = 0,1,...}, {uk £ 
W, fc = 0,l,...},and{dfc G RP, k ^ 0,1, ...} that satisfy 



Zfe+i - argmin/i(z) + ^||Gz-Ufc-dfc||2 
z 2 

Ufc+i - argmin/2(u) + ^||GzA;+i-u-dfe||2 
u 2 



< 
< 



and 



ife+i 



= dfe - (Gzfe+i - Ufc+i). 



Pk 
(16) 



Then, if jTSl) has a solution, say z* , the sequence {z^} con- 
verges to z* . If (T^ does not have a solution, then at least one 
of the sequences {u^} or {dk} diverges. 



According to Theorem[T] it is not necessary to exactly solve 
the minimizations in lines 3 and 4 of ADMM: as long as 
the sequences of errors are absolutely summable, convergence 
is not compromised. As shown in Section IIV-DI this fact is 
quite relevant in designing instances of ADMM, when these 
minimizations lack closed form solutions. 

The proof of Theorem [T] is based on the equivalence 
between ADMM and the DRS method applied to the dual 
of problem ( fTSb . For recent and comprehensive reviews of 
ADMM, DRS, and their relationship with Bregman and split- 
Bregman methods, see ifTsl . 061 . 

B. A Variant of ADMM 

Notice that the ADMM and the associated convergence 
theorem presented in the previous subsection apply to objec- 
tive functions of the form (flSl l. i.e., which are the sum of 
two functions. The fact that our objective functions (|9]l, ( fTTT i. 
and ( fTST l involve more than two terms raises the following 
question: what is the best way of mapping an objective with 
more than two terms into (flSl l so that the resulting ADMM is 
easily applicable and the conditions of Theorem [T] still hold. 
In this section, we give an answer to this question, which will 
constitute the core of our approach. 

Consider a generalization of problem (flSl l. where instead of 
two functions, we have J functions, that is. 



min V5,(HWz), 



(17) 



where gj : R'°j ^> R are closed, proper, convex functions, 
and H^-'^ G RPj^'' are arbitrary matrices. The minimization 
problem ( fTTI l can be written as ( fTsT i using the following 
correspondences: /i — 0, 

H(i) 



G 



G 



(18) 



where p = pi 



^pxd 



■ + pj, and /2 
/2(u) = ^g,(u(^')), 



given by 



(19) 



where u^^') G R^^' and u = [(u(i))'^, . . . , (u^'^))'^]'^ G R^. 

We are now in position to apply ADMM. The resulting 
algorithm has exactly the same structure as the one in Fig. [T] 
with 











dfc 









The fact that /i = turns Step 3 of the algorithm into a 
simple quadratic minimization problem, which has a unique 
solution if G has full column rank: 

|2 



argmin G z — Cfclh 



(G^G)-^G-Cfe 



j=l 



-1 J 



(20) 

^k ' 
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where (^j. = Uk + dk (and, naturally, C 



d[^'^) and 



the second equality results from the particular structure of G 
in ([H. 

Furthermore, our particular way of mapping problem (fTTI i 
into problem ( fTSl l allows decoupling the minimization in Step 
4 of ADMM into a set of J independent ones. In fact, 

Uk+i ^ argmin /2(u) + ^\\Gzk+i u dk\\l (21) 
u 2 

which can be written as 



u 



(1) 

k+l 



^k+1 



arg min gi{u^'>) + ■ ■ ■ + gj{u'^-^>) 

























Zfe+i - 






















. 4'^ _ 





Clearly, the minimizations with respect to u^^\ . . . , u*^''^ are 
decoupled, thus can be solved separately, leading to 



(7) • / \ A* II (l) l|2 

,yj> ^arg mm 5j(v) + - v-s^^ 



"fe+1 

for j = 1 , . . . , J, where 



(22) 



Eq. ( l22b defines the so-called Moreau proximity operators ||9] 
of gi, ■■■,gj, applied to s'j}\ s^,'''', respectively, denoted as 



,0-) 



(23) 



Some comments on the algorithm are relevant. Firstly, 
being exactly an ADMM, and since all the functions gj, for 
j = 1,...,J, are closed, proper, and convex, convergence is 
guaranteed if G has full column rank. This full column rank 
condition, which is also required for the inverse in (|20] | to 
exist, will be studied in the next section for each of the specific 
problems considered in this paper. 

For some functions, this mapping can be computed exactly 
in closed form. For example, if gj(x) = ||x||i, the correspond- 
ing proximity operator g-/^ is simply a soft threshold, 

*g,/M(v) = soft(v, = sign(v) max{|v| - (1/^), 0}, 

(24) 

where sign(-) denotes the component-wise application of the 
sign function, denotes the component- wise product, |v| 
denotes the vector of absolute values of the elements of v, and 
the maximum is computed in a component-wise fashion. For 
other functions, the corresponding Moreau proximity operator 
does not have a simple close form solution and needs to be 
computed numerically. 

IV. PoissoNiAN Image Reconstruction with 
TV-Based Regularization 

A. Applying ADMM 

In this section, we apply the algorithmic framework pre- 
sented in Section IIII-BI to the total-variation-based criterion 



Algorithm Poisson Image Decomolution by AL (PIDAL-TV) 

1. Choose u'^^ u'^\ d'^\ dj;^', d'^\ ^j, and r. Set fc <- 0. 


2. 

a 


repeat 


' "fe + °fc 




4. 








5. 








6. 


Ik ^ 






7. 




-(K^K + 2I)-S, 




8. 




-Kzfc+i-d^'' 


m 


9. 




<- argmin ^llv - i^l^'Hi + 




10. 


k 


- z? 1 1 — d. 


i = l 


11. 


(2) 

11 

"fc+1 


/i (2) ii2 
^ cLij^iiiiii 2 H k II 


rTVfv") 
; 1 V i^v;. 


12. 


'^k ^ 


j(3) 
— 7,1 11 — ri 




13. 


(3) 
"ic+l 


4— arg min — llv — u^J'^ |P + 
V 2 




14. 


"fe+1 


<-d«-(Kz,+i-uWj 




15. 


"fe+1 


j(2) / (2) X 




16. 




^df -(Zfe+i-u(^,) 




17. 


k -t^ k 


+ 1 




18. 


until some 


Stopping criterion is satisfied. 





Fig. 2. The PIDAL-TV algorithm. 

(l9]l-(fT0]i. The objective function in ( fTOl i has the form ( fTTI ) with 

J = 3, 

gi=jC-, 52 = T TV, ee ta:;: (25) 

and 

H(i) = K, H(2) = I, H^'^) = I. (26) 

The resulting ADMM algorithm, which we call PIDAL-TV 
(Poisson image deconvolution by augmented Lagrangian - 
total variation), is shown in Fig. |2l 



B. Implementation Aspects and Computational Cost of 
PIDAL-TV 

Notice that line 7 of PIDAL-TV corresponds to for the 
particular form of matrix G in this problem: G = [K^ I I]^ 
(see (|25] | and (|26]|). which is of course of full column rank. 
Moreover, if K models a periodic convolution, it is a block 
circulant matrix and the inversion in line 7 of the algorithm 
can be implemented in O(nlogn) operations, via the FFT 
algorithm. Although this is a well-known fact, we include the 
derivation in the next paragraph, for the sake of completeness. 

Assuming that the convolution is periodic, K is block- 
circulant with circulant blocks and can be factorized as 



K 



U^DU, 



(27) 



where U is the matrix that represents the 2D discrete Fourier 
transform (DFT), = U^^ is its inverse (U is unitary, i.e., 
UU^ = U^U = I), and D is a diagonal matrix containing 
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the DFT coefficients of the convolution operator represented 
by K. Thus (with = K^, since K is a real matrix) 

+ = (U^D*DU + 2U'^U)"^ (28) 

= (|D|2 + 2I)"^U, (29) 

where (•)* denotes complex conjugate and |Dp the squared 
absolute values of the entries of D. Since |Dp+2 1 is diagonal, 
its inversion has 0{n) cost. Products by U and can be 
carried out with O(nlogn) cost using the FFT algorithm. 

The minimization in line 9 of PIDAL-TV is separable. With 
respect to each component, it has the form 

mm|w + LR_^{v) -y log(u+) + ^{v ~ vf^ . (30) 

It is simple to show that the solution of ( l30l ) leads to 



'■i.k+l 




(31) 



where v^^l denotes the i-th component of ^'[.^^ Notice that 
,(1) 



'■i k+i is always a non-negative quantity. 

The minimization in line 1 1 of PIDAL-TV is, by definition, 
the Moreau proximity operator \I>(7./^)jv : R" ||9l, 

("2") 

which corresponds to applying TV-denoising to v^, '. Below, 
we address in detail the issue of how to implement this 
operator and its implications to the convergence of PIDAL- 
TV. Suffice it to say here that most TV-denoising algorithms 
have 0(ri) cost. 

The minimization in line 13 of PIDAL-TV corresponds to 

(3) 

the projection of v\ onto the first orthant, thus 



(32) 



where the maximum is to be understood in a component-wise 
sense; this projection has of course 0(n) cost. 

From the observations in the previous paragraphs, the com- 
putational costs of the lines of PIDAL-TV are the following. 
Lines 3, 4, 5, 9, 10, 11, 12, 13, 15, and 16 have 0{n) 
cost. Lines 6, 7, 8, and 14 have 0(n log n) cost. Thus the 
computational cost of PIDAL-TV scales as O(nlogn). 

C. Convergence of PIDAL-TV: Exact TV 

Convergence of PIDAL-TV is addressed by the following 
corollary of Theorem [T] for the (ideal) case where ^(t-/^)jv 
(line 11) is computed exactly. The minimizations in lines 9 
and 13 have the closed-form solutions given in ( l3Tl i and ( |32| |. 

Corollary 1: If the minimizations in lines 9, 11, and 13 of 
PIDAL-TV are solved exactly, then the algorithm converges to 
a minimizer of ifTOl). provided one exists. 

Proof: PIDAL-TV is an instance of ADMM in Fig. [T] 
where /i = and /2 has the form ( fT9l ). with J = 3 and the 
gi functions given in ( |25] |. which are all closed, proper, and 
convex. Function /2 is thus also closed, proper, and convex. 
Matrix G — [K."^ 11]^ has full column rank. The minimiza- 
tion in line 4 of ADMM corresponds to lines 9, 11, and 13 
of PIDAL-TV; if these minimizations are solved exactly, then 



according to Theorem \T\ convergence to a minimizer of the 
objective function, if one exists, is guaranteed. ■ 



D. Convergence of PIDAL-TV: Approximate TV 

As is well known, the TV denoising problem has no 
closed form solution, and many iterative algorithms have been 
proposed to solve it (see |4], 16), 110|, |31|, and references 
therein). Here, we adopt Chambolle's algorithm [41. 

Of course, in practice, Chambolle's (or any other iterative) 
algorithm can only run for a finite number of iterations, thus 
the minimization in line 1 1 of PIDAL-TV can only be solved 
approximately. However, as stated in Theorem [H this will 
not compromise the convergence of ADMM/PIDAL-TV, if 
the corresponding error sequence is summable. To achieve 
this goal, we adopt a simple procedure in which the internal 
variables of Chambolle's algorithm (the discrete gradient, see 
|4|) are initialized, in each iteration of PIDAL-TV, with those 
obtained in the previous iteration. We will now formalize this 
idea and provide experimental evidence that this procedure 
does produce a summable error sequence. 

Let us define /? = t/h and let (s, q) — ^^^{r, p) be the 
result of running t iterations of Chambolle's algorithm with 
its internal variables initialized at p, where s is the obtained 
(denoised) image (which is approximately 4'^Tv(r)) and q 
the final values of the internal variables. Consider now two 
possible implementations of line 1 1 of PIDAL-TV: 



11(a). 
11(b). 



, (2. a) N 



,T,(*) / (2) N 



(33) 



*?TV(4'\P)- (34) 



Implementation 11(a) uses the proposed internal variables ini- 
tialization, whereas in 11(b) the internal variables are always 
initialized to the same values (usually zeros). Consider now 
the corresponding error sequences 



u 



(2.a) 
fe+1 



*^Tv(«-i'^) 



„(2:b) ^ . (2)^ 



(35) 
(36) 



Notice that since the two other minimizations (lines 9 and 13) 
are solved exactly, the sequences p^,"-* and pf^ correspond to 
the sequence pk in Theorem [T] 

The following experiment provides evidence that p^"'' is 
summable, but p^^f^ is not. Consider the same setup as in the 
first experiment in fyf\: the original image is a portion of 
the Cameraman image scaled to a maximum value of 3000 
and then blurred with a Gaussian kernel of unit variance; the 
observed image is generated according to ([Til. As in |37|, we 
set T = 0.008 and p — r/50. The number of iterations 
of Chambolle's algorithm is set to 5 or 20. To compute 
'4'/3Tv(^'fe ) (almost) exactly, we run 4000 iterations of Cham- 
bolle's algorithm. In Figure[3] it is clear that the p^^^ sequences 
are not even decreasing, let alone summable. In contrast, the 
sequences p'"^^ approach zero, for both choices of t. Evidence 
for the summability of the p'jf' sequences is provided by the 
fact that by fitting a function of the form A{l/k)'^ to the 
tails of these sequences (i.e., for k = 20, ...,200), we obtain 
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Iterations (k) 



Fig. 3. Error sequences p^^"' and p^*' for t = 5 and t = 20 (number of 
iterations of Chambolle's algoritlim) and fitted functions of tlie form A to 
the sequences {p'j^Kk = 20, 21, ...200}. 



values of cu that are larger than one (cu ~ 1.3, for t — 5, and 
uj ~ 1.43, for t = 20). 

In conclusion, the experiment reported in the previous 
paragraph, though of course not a formal proof, strongly 
suggests that by implementing line 1 1 of PIDAL-TV as in ( [33] ), 
the corresponding error sequence (with respect to the exact 
minimizations) is summable, thus we can invoke Theorem [T] 
to state that PIDAL-TV converges. Moreover, this experiment 
shows that this is achieved with a quite small number of 
iterations in each call of Chambolle's algorithm. In all our 
experiments with PIDAL-TV, we thus use ( [33] l, with t = 5. 

V. PoissoNiAN Image Reconstruction with 
Frame-Based Regularization 

We now consider the frame-based analysis criterion ( fTTT ). 
and the frame-based synthesis criterion ( fT3] l. 

A. Analysis Criterion 

In this case, the objective function is given by ( fTSl l. which 
has the form ( fTTI l with J = 3, 



and 



91= -C, g2=T||-||i, 53 = 



(37) 



The resulting instance of ADMM, which we call PIDAL-FA 
(where FA stands for "frame analysis"), is shown in Fig. |4] 
The matrix being inverted in line 7 results from assuming that 
P is the analysis operator of a 1 -tight (Parseval) frame, thus 
(H(2))Th(2) = p^p = I. Notice that Hne 7 of PIDAL-FA 
corresponds to (|20] | for the particular form of matrix G in this 
case, G = [K^ P^ I]^, which of course has full column 
rank. As in PIDAL-TV, if K models a convolution, the inverse 
+ can be computed with 0{n\ogn) cost, using 

the FFT (see 

For most tight frames used in image processing, products 
by P and P^ correspond to the inverse and direct transforms 



Algorithm Poisson Image Deconvolution by AL (PIDAL-FA) 



I. 


Choose ul''', ul^\ 




2. 


repeat 




3. 






4. 








5. 








6. 


Ik ^ 




7. 




f- (K^ 


K + 2 


8. 




-Kzfc 





A3) 



(1) 

fc 



(2) 



1 

' 2 



k + 1 



(2) 



2,k 



■ /^ii (2)ii2 , n II 

arg mm - V - i^' ' +r v i. 

V 2 



j(3) 



j(2) 
.(3) 



j(2) 
.(3) 



Pz 



fc+1 



fc + 1 

^(2) 
-^k + 1 



C+1 ^ - ^^+1 

fc := A; + 1 



18. until some stopping criterion is satisfied. 



Fig. 4. The PIDAL-FA algorithm. 



for which fast algorithms exist. For example, in the case of 
translation-invariant wavelet transforms, these products can 
be computed using the undecimated wavelet transform with 
0{n\ogn) cost |25J, [28J. Curvelets also constitute a Parseval 
frame for which fast 0{n\ogn) implementations of the for- 
ward and inverse transform exist f3l. Yet another example of 
a redundant Parseval frame is provided by complex wavelets, 
with the corresponding direct and inverse transforms having 
0{n) cost [24], f34]. In conclusion, for a large class of 
choices of P, the cost of lines 6, 10, and 15 of PIDAL-FA is 
0{n\ogn). 

The expressions in lines 9 and 13 of PIDAL-FA are similar 
to those in lines 9 and 13 of PIDAL-TV, respectively; see also 
(©, (EB, and m. 

The minimization in line 11 is, by definition, the Moreau 
proximity operator of the £i norm [9J, which corresponds to 
a soft-threshold (l24l l. 

In summary, from the observations in the previous para- 
graphs, the computational costs of the lines of PIDAL-FA are 
the following. Lines 3, 4, 5, 9, 11, 12, 13, and 16 have 0{n) 
cost. Lines 6, 7, 8, 10, 14, and 15 have 0{n\ogn) cost. Thus 
the computational cost of PIDAL-FA scales as O(nlogn). 

Finally, convergence of PIDAL-FA is addressed by the 
following corollary of Theorem [T] 

Corollary 2: The PIDAL-FA algorithm converges to a min- 
imizer of ifTTI). provided one exists. 

Proof: The proof is similar to, but simpler than, that of 



Algorithm Pais son Image Deconvolution by AL (PIDAL-FS) 



1 . Choose u, 



(1) „(2) „(3) .(1) .(2) ,(3) 







d(\ dl^ ',d};'', n, and t. Set k ■ 



2. 


repeat 




a 
J. 


^fc ^ "fe +°k 




4. 


^(2) ^ (2) . .(2) 




5. 






6. 






7. 


Zfc+i (W^K^K W + I 


+ W^W) 


8. 







(2) 

"fc+1 
"fc+1 

d(2) 

d(3) 
"fe+1 



1 

" 2 



q{2) 



:{'^i^0}. 



10. 

11. 

12. 

13. 

14. 

15. 

16. 
17. 

18. until some stopping criterion is satisfied. 



d 

d^ 

■■k + l 



(2) 
fc 

(3) 



^(2) 



Fig. 5. The PIDAL-FS algorithm. 



Corollary [T] since all the minimizations involved are solved 
exactly in closed form. Clearly, matrix G = [K^ P-'^ I]-'^ 
has full column rank, thus Theorem [T] guarantees convergence 
to a minimizer of the objective function. ■ 



B. Synthesis Criterion 

In the synthesis formulation, the objective function is given 
by (fT3] l, which has the form ( fTTI l with J = 3, 



and 



5i=£, g2=Tl|-||i, 53 



The resulting ADMM algorithm, which we call PIDAL-FS 
(where FS stands for "frame synthesis"), is shown in Fig. |5] 
Notice that line 7 of PIDAL-FS corresponds to for 
the particular form of matrix G in this problem: G = 
[(KW)-'" I W-^]-^. This matrix has of course full column 
rank. However, even if K models a periodic convolution (thus 
is block circulant), the question remains of how to efficiently 
compute the matrix inverse in line 7, since K W is not block 
circulant. The next paragraph shows how to sidestep this 
difficulty. 

Consider that matrix W corresponds to a 1 -tight (Parse- 
val) frame, i.e., WW^ = I, and start by noticing that 
W^K^KW + I + W^W = W^(K^K + I)W + L 
Applying the Sherman-Morrison-Woodbury (SMW) matrix 



TABLE I 

Initialization of the PIDAL algorithms. 







(2) 










PIDAL-TV 


y 


y 


y 











PIDAL-FA 


y 


Py 


y 











PIDAL-FS 


y 




K'y 












inversion formula yields 



W^(K^K + I)W + I) ' = 



= I- W^fwW^+ (K^K + I) ^ 



W 



I W^(I+ (K"K + I 



W. 



(38) 



Using the factorization dZTb . we have 

I+(K^K + I)"')"' =U^(I+(|D|2+I)-i)"'u, (39) 

where both inversions have 0{n) cost since |Dp and I 
are diagonal, thus products by the matrix in ( [39] l have the 
0(nlog?i) cost associated to the FFT implementation of the 
products by U and U^. 

The leading cost of line 7 of PIDAL-FS (given by ( |38] l) 
will thus be either 0{n\ogn) or the cost of the products 
by and W. As mentioned above, for a large class of 

choices of frames, matrix-vector products by W and 
have 0(71 log n) cost. 

From the observations in the previous paragraphs, the com- 
putational costs of the lines of PIDAL-FS are the following. 
Lines 3, 4, 5, 9, 10, 11, 12, 13, 15, and 16 have 0[n) 
cost. Lines 6, 7, 8, and 14 have O(nlogn) cost. Thus the 
computational cost of PIDAL-FS scales as 0{n\ogn). 

Finally, convergence of PIDAL-FS is addressed by the 
following corollary of Theorem [T] 

Corollary 3: The PIDAL-FS algorithm converges to a mini- 
mizer of iTOl). provided one exists. 

Proof: The proof is similar to that of Corollary|2] since all 
the minimizations involved are solved exactly in closed form. 
Clearly, matrix G = [W^K^ I W^]^ has full column 
rank, thus Theorem [T] guarantees convergence to a minimizer 
of the objective function. ■ 



VI. Experiments 

We now report experiments where PIDAL is compared with 
other state-of-the-art methods, namely those proposed in 1 12 1, 
Iil9j , |i37|. All the algorithms are implemented in MATLAB 
and the experiments are carried out on a PC with a 3.0GHz 
Intel Core2Extreme CPU, with 4Gb of RAM, under Microsoft 
Windows Vista. Unless otherwise indicated, we adjust the reg- 
ularization parameter r to achieve the highest improvement in 
signal-noise-ratio (ISNR = 10 log^Q (||y — xjH/llx — x|||)). 
The PIDAL algorithms are initialized as shown in Table |T] 

According to Theorem [T] ADMM (thus PIDAL) converges 
for any choice of ji > {). However, this parameter does 
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influence the speed of the algorithms. To our knowledge, there 
is work on methods to choose this parameter for optimal speed; 
in our experiments, we use the following rule of thumb, found 
to achieve satisfactory results: /i = 60r/Af, where M is the 
maximum intensity of the original image. We have observed 
that the results do not change significantly if this parameter is 
changed to one order of magnitude below or above this choice. 



A. Comparison with l[37]l 

We begin by comparing with the algorithms (PIDSplit and 
PIDSplit+) proposed in ll37l . which (as acknowledged by the 
authors of Ii37j ) is based on the earlier version of PIDAL- 
TV ifm . The setup was already described in Section BV-DI 
the original image is a portion (84 x 84) of the Cameraman 
image, scaled to a maximum value of 3000 and blurred with 
a Gaussian kernel of unit variance; the observed image is 
generated according to ([T]l; the regularization parameter is 
set to T = 0.008. In the experiments reported in 1,37] , the 
TV denoising step of PIDAL-TV is implemented by an inner 
iterative algorithm with a tight stopping criterion based on the 
change between two consecutive images. Our implementation 
of PIDAL-TV, as explained in Section IIV-DI uses a small and 
fixed number of iterations (just 5) of Chambolle's algorithm, 
which is initialized as explained in that section. Because PID- 
Split and PIDSplitH- have no inner loop, each of its iterations 
is roughly equivalent to one iteration of PIDAL-TV with just 
one iteration of Chambolle's algorithm. In |37J, PIDSplit and 
PIDSplitH- were run for 2150 iterations; we thus run PIDAL- 
TV for 2150/5 = 430 iterations, corresponding to roughly the 
same amount of computation. Fig. |6]shows the evolution of the 
mean absolute error (MAE = ||x— x||i/n) and ISNR along the 
first 160 iterations of PIDAL-TV (as well as elapsed time); it is 
clear that convergence is achieved after less than 140 iterations 
(4.3 seconds, in our computer). This is dramatically less than 
what is reported in |37| for PIDAL-TV; in terms of iterations 
of PIDSplit and PIDSplitH-, this corresponds to approximately 
150 * 5 = 750 iterations, thus also much less than the 2150 
iterations (11 seconds) reported in that work. Finally, Fig. [T] 
shows the original, observed, and restored images; as expected, 
the image estimates produced by PIDSplit and PIDAL-TV are 
very similar. 

FinaUy, we also tested PIDAL-FA and PIDAL-FS on the 
same example, using a fully redundant Haar frame. The plots 
of ISNR and MAE are presented in Figs. [8] and |9] while the 
estimated images are shown in Fig. [TO] 

These results show that, in this example, PIDAL-FA per- 
forms slightly better than PIDAL-TV in terms of ISNR and 
similarly in terms of MAE, with PIDAL-FA achieving its 
best estimate faster than PIDAL-TV. The synthesis-based 
criterion implemented by PIDAL-FS is a little worse in terms 
of both ISNR and MAE, and PIDAL-FS also takes longer 
than PIDAL-FA to achieve its best estimate. This poorer 
performance of the synthesis formulation (in line with recent 
results in B51 ) was also found in all the experiments reported 
below, so we will only present results from PIDAL-TV and 
PIDAL-FA. 



time (seconds) 



, 0.5 1 l.S 2 2.5 3 3.5 4 4.5 




20 40 60 SO 100 120 140 160 
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Fig. 6. Evolution of tlie mean absolute error (MAE) and improvement in 
signal-noise-ratio (ISNR) along tlie iterations and elapsed time of PIDAL-TV, 
for the experiment of Section IVI-AI 




Fig. 7. Experiment of Section FVI-AI Top row: original (left); blun'ed and 
noisy image (right). Bottom row: estimate from Ii37] ; estimate by PIDAL-TV 
(ISNR=4.8dB). 
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Fig. 8. Experiment of Section IVl-AI Evolution of the mean absolute error 
(MAE) and improvement in signal-noise-ratio (ISNR) along the iterations and 
elapsed time of PIDAL-FA. 
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Fig. 9. Experiment of Section IVI-AI Evolution of the mean absolute error 
(MAE) and improvement in signal-noise-ratio (ISNR) along the iterations and 
elapsed time of PIDAL-FS. 




Fig. 10. Experiment of Section IvFaI Left: PIDAL-FA estimate (ISNR = 
5.3dB). Right: PIDAL-FS estimate (ISNR = 4.3dB). 

B. Comparison with /TP/ 

The next experiment follows f{9\: the original image is the 
complete (256 x 256) Cameraman, scaled to a maximum value 
of 17600, the blur is 9 x 9 uniform. As in the experiment 
reported in the previous subsection, this is a high SNR 
situation. Fig. [TT] shows the evolution of the MAE and ISNR 
along the execution of PIDAL-TV; it is clear that convergence 
is achieved after about 160 iterations (25 seconds, in our com- 
puter). A detail of the blurred, and estimated images (from 1 19] 
and using PIDAL-TV and PIDAL-FA) are shown in Fig. [12] 
Although the TV and FA regularizers are considerably simpler 
than the locally adaptive approximation techniques used in 
|fT9l . both PIDAL-TV and PIDAL-FA achieve higher ISNR 
values (7.0dB and 6.95dB, respectively) than that reported in 
lfr9J (6.61dB). 

C. Comparison with l[12]l 

In the last set of experiments we compare our approach with 
another recent state-of-the-art algorithm (herein referred to as 
DFS), proposed in lfT2l . for which the MATLAB implemen- 
tation is available at www.greyc.ensicaen.fr/~fdupe. 
That work includes comparisons with other methods, namely: 
Richardson-Lucy with multi-resolution support wavelet reg- 
ularization (RL-MRS) |[39l ; fast translation invariant tree- 
pruning reconstruction (FTITPR) li44J : Richardson-Lucy with 
total variation regularization (RL-TV) ifTTl . The results in llT2l 
show that the algorithm therein proposed generally achieves 
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Fig. 1 1 . Experiment of Section IVI-BI Evolution of the mean absolute error 
(MAE) and improvement in signal-noise-ratio (ISNR) along the iterations and 
elapsed time of PIDAL-TV. 




Fig. 12. Experiment of Section fVI-B I Top row: blurred noisy image (left) 
and estimate from |19| (ISNR=6.61dB). Bottom row: PIDAL-TV estimate 
(left, ISNR = V.OdB); PIDAL-FA estimate (right, lSNR=6.95dB). 

better performance (i.e., lower MAE) than the others, except 
for one of the images (a microscopy cell image) where RL- 
MRS outperforms DFS. For this reason, we will report results 
comparing PIDAL-TV and PIDAL-FA versus DFS and RL- 
MRS. For PIDAL-FA, we use a redundant Haar frame for the 
Cameraman image and Daubechies-4 for the other images. As 
in [12,1 . the original images are scaled to a maximum value 
M, belonging to {5, 30, 100, 255}, and then blurred by a 7 x 7 
uniform filter. 

The DFS algorithm does not include a stopping criterion, 
with the results reported in lfT2l having been obtained by 
running a fixed number (200) of iterations. In order to compare 
the running times of PIDAL-TV, PIDAL-FA, and DFS, we run 
DFS until the MAE decreases less than 0.01% between two 
consecutive iterations. Our algorithms are stopped when the 
following condition is met: 

|zfc - Zfc-l||2 ^ ^ 
l|Zfc-l||2 

with S = 0.005 if M = 5 and S = 0.001 in all the other cases. 
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Notice that this favors DFS, since a stopping criterion based 
on MAE is not appUcable in practice due to the absence of 
the original image. 

The results reported in Table show that, in 9 out of the 
12 experiments, either PIDAL-TV or PIDAL-FA achieves the 
lowest MAE. Notice however, that the main goal of this paper 
was not to introduce a new restoration criterion aiming at 
obtaining the lowest possible MAE (or any other performance 
measure), but rather to introduce algorithms to solve the 
optimization problems resulting from variational formulations 
of Poissonian image restoration. In terms of computational 
efficiency, PIDAL-TV and PIDAL-FA are clearly faster than 
the DRS algorithm, except in the very low SNR situations 
(M = 5) for two of the images (Cameraman and Cell). 

VII. Concluding Remarks 

We have propose new algorithms to handle the optimiza- 
tion problems resulting from regularization approaches to the 
restoration of Poissonian images. These optimization problems 
include several difficulties: the Poisson log-likelihood is non- 
quadratic and its gradient is not Lipschitz; the state-of-the- 
art regularizers are non-smooth; there is a non-negativity 
constraint. We have started by presenting sufficient conditions 
for existence and uniqueness of solutions of these optimization 
problems, for the following regularizers: total-variation, frame- 
based analysis, and frame-based synthesis. These problems 
were handled by adapting the alternating direction method 
of multipliers (ADMM) to their particular forms. This adap- 
tation is based on a new way of using ADMM to deal 
with problems in which the objective function is a linear 
combination of convex terms, which can be used in many other 
problems. We gave sufficient conditions for convergence and 
proved that these are met in the considered cases. Finally, 
we have experimentally compared the proposed algorithms 
against competing techniques, showing that they achieve state- 
of-the-art performance both in terms of speed and restoration 
accuracy. 

Appendix A: Convex Analysis 

We very briefly review some basic convex analysis results 
used in this paper. For more details see ||9l, ll47l . 

Consider a function / : A" — M = R U {— oo, +oo}, where 
M is called the extended real line, and A" is a real Hilbert 
space. The domain of function / is the set dom(/) — {x : 
/(x) < +oo}. 

The function / is convex if /(au + (1 — q;)v) < q;/(u) + 
(1 — q;)/(v), for any u, v e A" and any a e [0, 1]. Convexity 
is said to be strict if the inequality holds strictly (<) for any 
u, V S dom(/) and a e]0, 1[. 

The function is called proper if it is not equal to +oo 
everywhere and is never equal to — oo. 

The function / is lower semi-continuous (Isc) at v if 

lim inf f(x) > ffv), 
<5\.0 xeB(v,5)'^ ^ ) - J \ h 

where -B(v, i5) = {x : ||x — v|| < 5} is the (5-ball around v, 
and II • II is the norm in the Hilbert space X. A function is 
called Isc if it is Isc at every point of its domain. 



A function / is called coercive if it verifies 
lim||xj|^oo /(x) — +00. Proper, Isc, coercive functions 
play a key role in optimization via the following theorem [9J : 

Theorem 2: If / is a proper, Isc, coercive, convex function, 
then the set argminxeA' /(x) is nonempty. 

The next theorem concerns strictly convex functions. 

Theorem 3: If / is a strictly convex function, the set 
argmiiixgA' /(x) possesses at most one element. 
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